The subject matter of this application is related to the subject matter in a co-pending non-provisional application by the same inventor(s) as the instant application and filed on the same day as the instant application entitled, xe2x80x9cUsing Tiling to Improve Performance in a Sparse Symmetric Direct Matrix Solver,xe2x80x9dhaving serial number 09/316,173, and filing date May 21, 1999.
1. Field of the Invention
The present invention relates to computer systems for performing sparse matrix computations. More particularly, the present invention relates to a method and an apparatus that uses a hybrid approach to efficiently perform a column modification (CMOD) operation in solving a system of linear algebraic equations involving a sparse coefficient matrix.
2. Related Art
The solution of large sparse symmetric linear systems of equations constitutes the primary computational cost in numerous applications, such as finite-element design, linear programming, circuit simulation and semiconductor device modeling. Efficient solution of such systems has long been the subject of research, and considerable progress has been made in developing efficient algorithms to this end. A direct solution technique known as xe2x80x9cCholesky Factorizationxe2x80x9d is the most widely used approach to solve such a system. Under Cholesky factorization, the complete solution sequence requires many stages, including matrix reordering, symbolic factorization, numerical factorization and triangular solution. Of these stages, numerical factorization is typically the most computationally expensive.
One method of performing numerical factorization is based on a right-looking supernode-supernode method described in xe2x80x9cParallel Algorithms for Sparse Linear Systemsxe2x80x9d by Michael T. Heath, Esmond Ng and Barry W. Peyton, in xe2x80x9cParallel Algorithms for Matrix Computationsxe2x80x9d by Gallivan, et al. (Editors), SIAM (1994) (referred to as HNP). In a sparse matrix, a supernode is a set of contiguous columns that have essentially the same sparsity structure. Supernodes can be used to organize the numerical factorization stage around matrix-vector (supernode-column) and matrix-matrix (supernode-supernode) primitive operations leading to a substantial performance improvement arising from more efficient use of the processor caches and pipelining units.
The numerical factorization step involves two fundamental sub-tasks:
(1) cdiv(j): division of column j of factor by a scalar; and
(2) cmod(j,k): modification of column j by column k, k less than j.
These sub-tasks can be organized around supernodes. For example, cdiv(j) can be organized as an internal factorization/update of supernode j, and cmod(j,k) can be organized as a modification of supernode j by supernode k, k less than j.
Typically, the second sub-task is where the bulk of the computational cost is incurred. In order to increase computational efficiency, the cmod(j,k) operation can be divided into three steps:
(a) computation of the update and accumulation into a temporary array;
(b) carrying out the non-zero index matching between the first columns of the source and destination supernodes and computing relative indices; and
(c) scattering updates from the temporary vector into the target destination supernode.
By dividing the cmod(j,k) operation in this way, it is possible to apply techniques used in dense matrix operations in the step (a). Note that step (a) is where the dominant amount of time is spent. In the discussion that follows, we refer the step (a) as the xe2x80x9clocal dense cmod operationxe2x80x9d. The local dense cmod operation involves computing a trapezoidal-shaped dense update that can be represented as a combination of a dense rank-k update and a matrix multiplication.
Library routines can be used to speed up the cmod computation. These library routines are typically written in assembly language and are hand-tuned for a specific machine architecture. For example, on the current generation of UltraSparc-II-based machines, the Sun Performance Library (see http://www.sun.com/workshop/performance/wp-perflib/) provides SPARC assembly-language implementations of BLAS1, BLAS2 and BLAS3 routines. These hand-tuned assembly language implementations can yield performance close to the theoretical peak of the underlying processor.
For example, portions of the cmod operation can be efficiently performed by invoking the BLAS3 xe2x80x9cdgemmxe2x80x9d matrix multiplication code from Sun Performance library. Unfortunately, invoking the BLAS3 dgemm matrix multiplication code requires supernodes to be copied into and out of temporary storage because of incompatibility between data-structures used by a typical sparse solver and those expected by the dgemm API. This copying can add a significant overhead. Consequently, using the BLAS3 xe2x80x9cdgemmxe2x80x9d matrix multiplication code only makes sense for supernodes above a certain size. Otherwise, the performance gains from using the BLAS3 xe2x80x9cdgemmxe2x80x9d library code are cancelled out by the additional overhead involved in copying.
Hence, what is needed is a system that performs the cmod operation using library routines in cases where the performance gains from using the library routines exceed the computational overhead required to use the library routines.
Another difficulty in attaining high performance in numerical factorization is due to the fact that supernodes can have varying shapes, sizes, and sparsity patterns. These varying shapes, sizes and sparsity patterns can greatly influence computational performance. In order to optimize computational performance for the cmod operation, the supernodes of varying shapes and sizes must be divided into smaller sub-units for computation so as to balance computational operations with memory references in a way that is tuned for the particular machine architecture on which the computation is to be run.
One embodiment of the present invention provides a system for efficiently perform a cmod operation in solving a system of linear algebraic equations involving a sparse coefficient matrix. The system operates by identifying supernodes in the sparse matrix, wherein each supernode comprises a set of contiguous columns having a substantially similar pattern of non-zero elements. In solving the equation, the system performs a CMOD operation between a source supernode and a destination supernode. As part of this CMOD operation, the system determines dimensions of the source supernode and the destination supernode. If a result of a function on the dimensions is lower than a threshold value, the system performs the CMOD operation between the source supernode and the destination supernode using a kernel that is written in an architecture-independent high-level language. If the result of the function on the dimensions is greater than the threshold value, the system performs the CMOD operation using a kernel that is written in an architecture-specific assembly code that runs efficiently on a specific computer architecture. Using this architecture-specific assembly code may require first copying supernodes into temporary data structures that are tailored for the assembly code.
In one embodiment of the present invention, the result of the function on the dimensions is lower than the threshold value if the number of columns of the source supernode is less than a first threshold value, the number of columns of the destination supernode is less than a second threshold value, and the number of common contiguous elements in columns to be updated in the destination supernode is less than a third threshold value. These threshold values can be empirically pre-computed using representative problem instances.
Thus, the present invention automates the selection of the appropriate kernel to maximize the overall performance of the cmod operation independently of the problem structure and size.